
/*
 * Author
 *   David Blom, TU Delft. All rights reserved.
 */

#include "gtest/gtest.h"
#include "SDCTubeFlowFluidSolver.H"
#include "SDCTubeFlowLinearizedSolidSolver.H"
#include "GaussRadau.H"
#include "GaussLobatto.H"
#include "ESDIRK.H"
#include "AdaptiveTimeStepper.H"
#include "SDCFsiSolver.H"
#include "AndersonPostProcessing.H"
#include "RBFCoarsening.H"
#include "RelativeConvergenceMeasure.H"
#include "ResidualRelativeConvergenceMeasure.H"
#include "MinIterationConvergenceMeasure.H"
#include "AitkenPostProcessing.H"
#include "AbsoluteConvergenceMeasure.H"

class SDIRKFsiSolidSolverTest : public ::testing::Test
{
    protected:
        virtual void SetUp()
        {
            scalar r0 = 3.0e-3;
            scalar h = 3.0e-4;
            scalar L = 0.126;
            scalar rho_s = 1000;
            scalar E0 = 4.0e5;
            scalar G = 4.0e5;
            scalar nu = 0.5;

            scalar a0 = M_PI * r0 * r0;
            scalar u0 = 0.26;
            scalar p0 = 0;
            scalar dt = 1;
            int N = 5;
            scalar T = 1;
            scalar rho_f = 1060;
            scalar E = 490;
            scalar cmk = std::sqrt( E * h / (2 * rho_f * r0) );

            bool parallel = false;
            int extrapolation = 0;
            scalar tol = 1.0e-3;
            int maxIter = 20;
            scalar initialRelaxation = 1.0e-3;
            int maxUsedIterations = 50;
            int nbReuse = 0;

            scalar singularityLimit = 1.0e-13;
            int reuseInformationStartingFromTimeIndex = 0;
            bool scaling = false;
            bool updateJacobian = false;
            scalar beta = 0.5;
            int minIter = 5;

            std::shared_ptr<tubeflow::SDCTubeFlowFluidSolver> fluid( new tubeflow::SDCTubeFlowFluidSolver( a0, u0, p0, dt, cmk, N, L, T, rho_f ) );

            std::shared_ptr<fsi::BaseMultiLevelSolver> solid;
            solid = std::shared_ptr<fsi::BaseMultiLevelSolver>( new tubeflow::SDCTubeFlowLinearizedSolidSolver( N, nu, rho_s, h, L, dt, G, E0, r0, T ) );

            assert( solid );

            shared_ptr<RBFFunctionInterface> rbfFunction;
            shared_ptr<RBFInterpolation> rbfInterpolator;
            shared_ptr<RBFCoarsening> rbfInterpToCouplingMesh;
            shared_ptr<RBFCoarsening> rbfInterpToMesh;

            rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
            rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
            rbfInterpToCouplingMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

            rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
            rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
            rbfInterpToMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

            shared_ptr<MultiLevelSolver> fluidSolver( new MultiLevelSolver( fluid, fluid, rbfInterpToCouplingMesh, rbfInterpToMesh, 0, 0 ) );

            rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
            rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
            rbfInterpToCouplingMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

            rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
            rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
            rbfInterpToMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

            shared_ptr<MultiLevelSolver> solidSolver( new MultiLevelSolver( solid, fluid, rbfInterpToCouplingMesh, rbfInterpToMesh, 1, 0 ) );

            std::shared_ptr< std::list<std::shared_ptr<ConvergenceMeasure> > > convergenceMeasures;
            convergenceMeasures = std::shared_ptr<std::list<std::shared_ptr<ConvergenceMeasure> > >( new std::list<std::shared_ptr<ConvergenceMeasure> >() );

            convergenceMeasures->push_back( std::shared_ptr<ConvergenceMeasure>( new RelativeConvergenceMeasure( 0, true, tol ) ) );
            convergenceMeasures->push_back( std::shared_ptr<ConvergenceMeasure>( new MinIterationConvergenceMeasure( 0, false, minIter ) ) );

            shared_ptr<MultiLevelFsiSolver> fsi( new MultiLevelFsiSolver( fluidSolver, solidSolver, convergenceMeasures, parallel, extrapolation ) );

            shared_ptr<PostProcessing> postProcessing( new AndersonPostProcessing( fsi, maxIter, initialRelaxation, maxUsedIterations, nbReuse, singularityLimit, reuseInformationStartingFromTimeIndex, scaling, beta, updateJacobian ) );

            std::shared_ptr<sdc::SDCFsiSolverInterface> sdcFluidSolver = std::dynamic_pointer_cast<sdc::SDCFsiSolverInterface>( fluid );
            std::shared_ptr<sdc::SDCFsiSolverInterface> sdcSolidSolver = std::dynamic_pointer_cast<sdc::SDCFsiSolverInterface>( solid );

            assert( sdcFluidSolver );
            assert( sdcSolidSolver );

            std::shared_ptr<fsi::SDCFsiSolver> fsiSolver( new fsi::SDCFsiSolver( sdcFluidSolver, sdcSolidSolver, postProcessing, extrapolation ) );

            std::shared_ptr<sdc::AdaptiveTimeStepper> adaptiveTimeStepper( new sdc::AdaptiveTimeStepper( false ) );

            std::string method = "ESDIRK53PR";

            esdirk = std::shared_ptr<sdc::ESDIRK> ( new sdc::ESDIRK( fsiSolver, method, adaptiveTimeStepper ) );
        }

        virtual void TearDown()
        {
            esdirk.reset();
        }

        std::shared_ptr<sdc::ESDIRK> esdirk;
};

TEST_F( SDIRKFsiSolidSolverTest, object )
{
    ASSERT_TRUE( true );
}

TEST_F( SDIRKFsiSolidSolverTest, timeStep )
{
    esdirk->solveTimeStep( 0 );
}

TEST_F( SDIRKFsiSolidSolverTest, run )
{
    esdirk->run();
}

TEST( SDIRKFsiSolidTest, linear )
{
    scalar r0 = 3.0e-3;
    scalar h = 3.0e-4;
    scalar L = 0.126;
    scalar rho_s = 1000;
    scalar E0 = 4.0e5;
    scalar G = 4.0e5;
    scalar nu = 0.5;

    scalar a0 = M_PI * r0 * r0;
    scalar u0 = 0.26;
    scalar p0 = 0;
    int N = 5;
    scalar T = 1;
    scalar rho_f = 1060;
    scalar E = 490;
    scalar cmk = std::sqrt( E * h / (2 * rho_f * r0) );

    bool parallel = false;
    int extrapolation = 0;
    scalar tol = 1.0e-8;
    int maxIter = 50;
    scalar initialRelaxation = 1.0e-3;
    int maxUsedIterations = 50;
    int nbReuse = 0;
    scalar singularityLimit = 1.0e-13;
    int reuseInformationStartingFromTimeIndex = 0;
    bool scaling = false;
    scalar beta = 0.01;
    bool updateJacobian = false;
    int minIter = 5;

    int nbComputations = 3;

    std::deque<std::shared_ptr<tubeflow::SDCTubeFlowFluidSolver> > fluidSolvers;
    std::deque<std::shared_ptr<tubeflow::SDCTubeFlowLinearizedSolidSolver> > solidSolvers;
    std::deque<int> nbTimeStepsList;

    for ( int iComputation = 0; iComputation < nbComputations; iComputation++ )
    {
        int nbTimeSteps = 120 * std::pow( 2, iComputation );
        std::cout << "nbTimeSteps = " << nbTimeSteps << std::endl;
        scalar dt = T / nbTimeSteps;

        std::shared_ptr<tubeflow::SDCTubeFlowFluidSolver> fluid( new tubeflow::SDCTubeFlowFluidSolver( a0, u0, p0, dt, cmk, N, L, T, rho_f ) );
        std::shared_ptr<tubeflow::SDCTubeFlowLinearizedSolidSolver> solid( new tubeflow::SDCTubeFlowLinearizedSolidSolver( N, nu, rho_s, h, L, dt, G, E0, r0, T ) );

        shared_ptr<RBFFunctionInterface> rbfFunction;
        shared_ptr<RBFInterpolation> rbfInterpolator;
        shared_ptr<RBFCoarsening> rbfInterpToCouplingMesh;
        shared_ptr<RBFCoarsening> rbfInterpToMesh;

        rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
        rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
        rbfInterpToCouplingMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

        rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
        rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
        rbfInterpToMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

        shared_ptr<MultiLevelSolver> fluidSolver( new MultiLevelSolver( fluid, fluid, rbfInterpToCouplingMesh, rbfInterpToMesh, 0, 0 ) );

        rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
        rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
        rbfInterpToCouplingMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

        rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
        rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
        rbfInterpToMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

        shared_ptr<MultiLevelSolver> solidSolver( new MultiLevelSolver( solid, fluid, rbfInterpToCouplingMesh, rbfInterpToMesh, 1, 0 ) );

        std::shared_ptr< std::list<std::shared_ptr<ConvergenceMeasure> > > convergenceMeasures;
        convergenceMeasures = std::shared_ptr<std::list<std::shared_ptr<ConvergenceMeasure> > >( new std::list<std::shared_ptr<ConvergenceMeasure> >() );

        convergenceMeasures->push_back( std::shared_ptr<ConvergenceMeasure>( new RelativeConvergenceMeasure( 0, true, tol ) ) );
        convergenceMeasures->push_back( std::shared_ptr<ConvergenceMeasure>( new MinIterationConvergenceMeasure( 0, false, minIter ) ) );

        shared_ptr<MultiLevelFsiSolver> fsi( new MultiLevelFsiSolver( fluidSolver, solidSolver, convergenceMeasures, parallel, extrapolation ) );

        shared_ptr<PostProcessing> postProcessing( new AndersonPostProcessing( fsi, maxIter, initialRelaxation, maxUsedIterations, nbReuse, singularityLimit, reuseInformationStartingFromTimeIndex, scaling, beta, updateJacobian ) );

        std::shared_ptr<sdc::SDCFsiSolverInterface> sdcFluidSolver = std::dynamic_pointer_cast<sdc::SDCFsiSolverInterface>( fluid );
        std::shared_ptr<sdc::SDCFsiSolverInterface> sdcSolidSolver = std::dynamic_pointer_cast<sdc::SDCFsiSolverInterface>( solid );

        assert( sdcFluidSolver );
        assert( sdcSolidSolver );

        std::shared_ptr<fsi::SDCFsiSolver> fsiSolver( new fsi::SDCFsiSolver( sdcFluidSolver, sdcSolidSolver, postProcessing, extrapolation ) );

        std::shared_ptr<sdc::AdaptiveTimeStepper> adaptiveTimeStepper( new sdc::AdaptiveTimeStepper( false ) );

        std::string method = "ESDIRK53PR";

        std::shared_ptr<sdc::ESDIRK> esdirk( new sdc::ESDIRK( fsiSolver, method, adaptiveTimeStepper ) );

        esdirk->run();

        fluidSolvers.push_back( fluid );
        solidSolvers.push_back( solid );
        nbTimeStepsList.push_back( nbTimeSteps );
    }

    std::cout << "solid" << std::endl;

    for ( int i = 0; i < 2; i++ )
    {
        fsi::vector ref;

        if ( i == 0 )
            ref = solidSolvers.back()->r;
        else
            ref = solidSolvers.back()->u;

        std::deque<scalar> errors;

        for ( int iComputation = 0; iComputation < nbComputations - 1; iComputation++ )
        {
            fsi::vector data;

            if ( i == 0 )
                data = solidSolvers.at( iComputation )->r;
            else
                data = solidSolvers.at( iComputation )->u;

            scalar error = (ref - data).norm() / ref.norm();
            errors.push_back( error );
        }

        for ( int iComputation = 0; iComputation < nbComputations - 2; iComputation++ )
        {
            scalar order = ( std::log10( errors.at( iComputation ) ) - std::log10( errors.at( iComputation + 1 ) ) ) / ( std::log10( nbTimeStepsList.at( iComputation + 1 ) ) - std::log10( nbTimeStepsList.at( iComputation ) ) );
            std::cout << "order = " << order << std::endl;

            if ( i == 0 )
                ASSERT_NEAR( order, 3, 0.1 );
        }
    }

    std::cout << "fluid" << std::endl;

    for ( int i = 0; i < 3; i++ )
    {
        fsi::vector ref;

        if ( i == 0 )
            ref = fluidSolvers.back()->u;

        if ( i == 1 )
            ref = fluidSolvers.back()->a;

        if ( i == 2 )
            ref = fluidSolvers.back()->p;

        std::deque<scalar> errors;

        for ( int iComputation = 0; iComputation < nbComputations - 1; iComputation++ )
        {
            fsi::vector data;

            if ( i == 0 )
                data = fluidSolvers.at( iComputation )->u;

            if ( i == 1 )
                data = fluidSolvers.at( iComputation )->a;

            if ( i == 2 )
                data = fluidSolvers.at( iComputation )->p;

            scalar error = (ref - data).norm() / ref.norm();
            errors.push_back( error );
        }

        for ( int iComputation = 0; iComputation < nbComputations - 2; iComputation++ )
        {
            scalar order = ( std::log10( errors.at( iComputation ) ) - std::log10( errors.at( iComputation + 1 ) ) ) / ( std::log10( nbTimeStepsList.at( iComputation + 1 ) ) - std::log10( nbTimeStepsList.at( iComputation ) ) );
            std::cout << "order = " << order << std::endl;

            if ( i == 1 || i == 2 )
                ASSERT_NEAR( order, 3, 0.1 );
        }
    }
}

TEST( SDIRKFsiSolidTest, linearized )
{
    scalar r0 = 3.0e-3;
    scalar h = 3.0e-4;
    scalar L = 0.126;
    scalar rho_s = 1000;
    scalar E0 = 4.0e5;
    scalar G = 4.0e5;
    scalar nu = 0.5;

    scalar a0 = M_PI * r0 * r0;
    scalar u0 = 0.26;
    scalar p0 = 0;
    int N = 5;
    scalar T = 1;
    scalar rho_f = 1060;
    scalar E = 490;
    scalar cmk = std::sqrt( E * h / (2 * rho_f * r0) );

    bool parallel = false;
    int extrapolation = 0;
    scalar tol = 1.0e-8;
    int maxIter = 50;
    scalar initialRelaxation = 1.0e-3;
    int maxUsedIterations = 50;
    int nbReuse = 0;
    scalar singularityLimit = 1.0e-13;
    int reuseInformationStartingFromTimeIndex = 0;
    bool scaling = false;
    scalar beta = 0.01;
    bool updateJacobian = false;
    int minIter = 5;

    int nbComputations = 3;

    std::deque<std::shared_ptr<tubeflow::SDCTubeFlowFluidSolver> > fluidSolvers;
    std::deque<std::shared_ptr<tubeflow::SDCTubeFlowLinearizedSolidSolver> > solidSolvers;
    std::deque<int> nbTimeStepsList;

    for ( int iComputation = 0; iComputation < nbComputations; iComputation++ )
    {
        int nbTimeSteps = 120 * std::pow( 2, iComputation );
        std::cout << "nbTimeSteps = " << nbTimeSteps << std::endl;
        scalar dt = T / nbTimeSteps;

        std::shared_ptr<tubeflow::SDCTubeFlowFluidSolver> fluid( new tubeflow::SDCTubeFlowFluidSolver( a0, u0, p0, dt, cmk, N, L, T, rho_f ) );
        std::shared_ptr<tubeflow::SDCTubeFlowLinearizedSolidSolver> solid( new tubeflow::SDCTubeFlowLinearizedSolidSolver( N, nu, rho_s, h, L, dt, G, E0, r0, T ) );

        shared_ptr<RBFFunctionInterface> rbfFunction;
        shared_ptr<RBFInterpolation> rbfInterpolator;
        shared_ptr<RBFCoarsening> rbfInterpToCouplingMesh;
        shared_ptr<RBFCoarsening> rbfInterpToMesh;

        rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
        rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
        rbfInterpToCouplingMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

        rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
        rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
        rbfInterpToMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

        shared_ptr<MultiLevelSolver> fluidSolver( new MultiLevelSolver( fluid, fluid, rbfInterpToCouplingMesh, rbfInterpToMesh, 0, 0 ) );

        rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
        rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
        rbfInterpToCouplingMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

        rbfFunction = shared_ptr<RBFFunctionInterface>( new TPSFunction() );
        rbfInterpolator = shared_ptr<RBFInterpolation>( new RBFInterpolation( rbfFunction ) );
        rbfInterpToMesh = shared_ptr<RBFCoarsening> ( new RBFCoarsening( rbfInterpolator ) );

        shared_ptr<MultiLevelSolver> solidSolver( new MultiLevelSolver( solid, fluid, rbfInterpToCouplingMesh, rbfInterpToMesh, 1, 0 ) );

        std::shared_ptr< std::list<std::shared_ptr<ConvergenceMeasure> > > convergenceMeasures;
        convergenceMeasures = std::shared_ptr<std::list<std::shared_ptr<ConvergenceMeasure> > >( new std::list<std::shared_ptr<ConvergenceMeasure> >() );

        convergenceMeasures->push_back( std::shared_ptr<ConvergenceMeasure>( new RelativeConvergenceMeasure( 0, true, tol ) ) );
        convergenceMeasures->push_back( std::shared_ptr<ConvergenceMeasure>( new MinIterationConvergenceMeasure( 0, false, minIter ) ) );

        shared_ptr<MultiLevelFsiSolver> fsi( new MultiLevelFsiSolver( fluidSolver, solidSolver, convergenceMeasures, parallel, extrapolation ) );

        shared_ptr<PostProcessing> postProcessing( new AndersonPostProcessing( fsi, maxIter, initialRelaxation, maxUsedIterations, nbReuse, singularityLimit, reuseInformationStartingFromTimeIndex, scaling, beta, updateJacobian ) );

        std::shared_ptr<sdc::SDCFsiSolverInterface> sdcFluidSolver = std::dynamic_pointer_cast<sdc::SDCFsiSolverInterface>( fluid );
        std::shared_ptr<sdc::SDCFsiSolverInterface> sdcSolidSolver = std::dynamic_pointer_cast<sdc::SDCFsiSolverInterface>( solid );

        assert( sdcFluidSolver );
        assert( sdcSolidSolver );

        std::shared_ptr<fsi::SDCFsiSolver> fsiSolver( new fsi::SDCFsiSolver( sdcFluidSolver, sdcSolidSolver, postProcessing, extrapolation ) );

        std::shared_ptr<sdc::AdaptiveTimeStepper> adaptiveTimeStepper( new sdc::AdaptiveTimeStepper( false ) );

        std::string method = "ESDIRK53PR";

        std::shared_ptr<sdc::ESDIRK> esdirk( new sdc::ESDIRK( fsiSolver, method, adaptiveTimeStepper ) );

        esdirk->run();

        fluidSolvers.push_back( fluid );
        solidSolvers.push_back( solid );
        nbTimeStepsList.push_back( nbTimeSteps );
    }

    std::cout << "solid" << std::endl;

    for ( int i = 0; i < 2; i++ )
    {
        fsi::vector ref;

        if ( i == 0 )
            ref = solidSolvers.back()->r;
        else
            ref = solidSolvers.back()->u;

        std::deque<scalar> errors;

        for ( int iComputation = 0; iComputation < nbComputations - 1; iComputation++ )
        {
            fsi::vector data;

            if ( i == 0 )
                data = solidSolvers.at( iComputation )->r;
            else
                data = solidSolvers.at( iComputation )->u;

            scalar error = (ref - data).norm() / ref.norm();
            errors.push_back( error );
        }

        for ( int iComputation = 0; iComputation < nbComputations - 2; iComputation++ )
        {
            scalar order = ( std::log10( errors.at( iComputation ) ) - std::log10( errors.at( iComputation + 1 ) ) ) / ( std::log10( nbTimeStepsList.at( iComputation + 1 ) ) - std::log10( nbTimeStepsList.at( iComputation ) ) );
            std::cout << "order = " << order << std::endl;

            if ( i == 0 )
                ASSERT_NEAR( order, 3, 0.1 );
        }
    }

    std::cout << "fluid" << std::endl;

    for ( int i = 0; i < 3; i++ )
    {
        fsi::vector ref;

        if ( i == 0 )
            ref = fluidSolvers.back()->u;

        if ( i == 1 )
            ref = fluidSolvers.back()->a;

        if ( i == 2 )
            ref = fluidSolvers.back()->p;

        std::deque<scalar> errors;

        for ( int iComputation = 0; iComputation < nbComputations - 1; iComputation++ )
        {
            fsi::vector data;

            if ( i == 0 )
                data = fluidSolvers.at( iComputation )->u;

            if ( i == 1 )
                data = fluidSolvers.at( iComputation )->a;

            if ( i == 2 )
                data = fluidSolvers.at( iComputation )->p;

            scalar error = (ref - data).norm() / ref.norm();
            errors.push_back( error );
        }

        for ( int iComputation = 0; iComputation < nbComputations - 2; iComputation++ )
        {
            scalar order = ( std::log10( errors.at( iComputation ) ) - std::log10( errors.at( iComputation + 1 ) ) ) / ( std::log10( nbTimeStepsList.at( iComputation + 1 ) ) - std::log10( nbTimeStepsList.at( iComputation ) ) );
            std::cout << "order = " << order << std::endl;

            if ( i == 1 || i == 2 )
                ASSERT_NEAR( order, 3, 0.1 );
        }
    }
}
